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ABSTRACT 

We extend our previous model-independent determination of the intergalactic background light, based purely on 
galaxy survey data, out to a wavelength of 5 /im. Our approach enables us to constrain the range of photon densities, 
based on the uncertainties from observationally determined luminosity densities and colors. We further determine 
a 68% confidence upper and lower limit on the opacity of the universe to y-rays up to energies of 1 .6/(1 + z ) TeV. 
A comparison of our lower limit redshift-dependent opacity curves to the opacity limits derived from the results of 
both ground-based air Cerenkov telescope and Fermi- LAT observations of PKS 1424+240 allows us to place a new 
upper limit on the redshift of this source, independent of IBL modeling. 
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1. INTRODUCTION 

Past work on estimating the spectral and redshift charac- 
teristics of the intergalactic photon density (IBL; with the 
z = 0 IBL usually referred to as the EBL) has depended 
on various assumptions as to the evolution of stellar pop- 
ulations and dust absorption in galaxies. There have also 
been attempts to probe the extragalactic background light 
(EBL) using studies of blazar y-ray spectra (Ackermann 
et al. 2012; Abramowski et al. 2013), an approach originally 
suggested by Stecker et al. (1992). 

We previously have pursued an independent, fully empiri- 
cal approach to calculating the IBL using deep galaxy survey 
data (Stecker et al. 2012 — hereafter SMS; see also Helgason & 
Kashlinsky 2012). In this paper, we extend the previous SMS 
results on the y-ray opacity of the universe to TeV energies 
accessible to ground-based air Cerenkov telescopes, here in- 
cluding galaxy survey data extending beyond the photometric 
I band to the M band (i.e., 0.8 /i m to ~5 /im). We stress that 
the SMS approach, being completely observationally based, 
is model independent, relying only on published luminosity 
densities. Such an approach has a particular advantage over 
model-based methods, as it enables a determination of both the 
IBL and its observational uncertainties without making any 
assumptions about how the galaxy luminosity functions (LFs) 
evolve, as some approaches require (e.g. Malkan & Stecker 
1998, 2001; Kneiske et al. 2002; Stecker et al. 2006; 
Franceschini et al. 2008; Finke et al. 2010; Kneiske & Dole 
2010) or by using semi-analytic models that make assumptions 
concerning galaxy evolution, stellar population synthesis mod- 
els, star formation rates, and dust attenuation (e.g., Gilmore et al. 
2009; Somerville et al. 2012). Dominguez et al. (201 1) used the 
redshift evolution of the K-band LFs of Cirasuolo et al. (2010), 
together with model templates based on AEGIS data to model 
the IBL. A recent detailed review of the problem has been given 
by Dwek & Krennrich (2013). 

Direct integration of the galaxy survey data in the wave- 
length range from the far-ultraviolet (FUV) to ~5 pm has now 
become possible. Observational data in this wavelength range. 


even out to redshifts 2 < z < 4, have become sufficiently 
complete to eliminate the necessity of using model-based ap- 
proaches. The high-redshift galaxies are sufficiently sampled 
to enable us to interpolate between observationally determined 
luminosity densities obtained for specific wavebands. We use 
observationally determined galaxy colors and their associated 
errors to fill in the gaps in wavelength bands where there are 
less data. 

Here we apply the technique of SMS to determine the IBL 
to a wavelength of ~5 /im, thus extending our previous range 
of y- ray opacities of the universe to energies up to 1.6 TeV. 
This range is not only relevant for the sources with the highest 
redshifts observed by the Fermi- LAT Telescope, but also for 
results in the sub-TeV and TeV energy range made using ground- 
based air Cerenkov telescopes. Our technique to determine the 
expected y-ray opacity of the universe is thus complementary to 
the technique of using y-ray observations directly, because the 
intrinsic (unabsorbed) spectra of the y -ray sources are uncertain. 

To further show the utility of our approach, we apply our 
lower limit opacities to the results of both ground-based air 
Cerenkov telescope observations and Fermi-L AT observations 
of the BL Lac object PKS 1424+240, which allows us to place 
a true upper limit on the redshift of this source, independent of 
IBL models. Previous observations of the ultraviolet absorption 
from the intergalactic medium of this object place a lower limit 
on its redshift of z ^ .6035 (Furniss et al. 2013). 

Our results give the y-ray opacity as a function of energy 
and redshift to within an observationally determined 68% 
confidence band. Therefore, a direct comparison with y-ray 
spectral data will allow an independent determination of the 
previously suggested effects on such spectra. These effects 
include secondary y-ray production generated by (1) cosmic- 
ray interactions along the line of sight to the source (Essey 
et al. 2010; Essey & Kusenko 2012) and (2) line-of-sight 
photon-axion oscillations during propagation (e.g., De Angelis 
et al. 2007). Both of these processes can produce an apparent 
reduction of the pair-production opacity effect as derived from 
y-ray spectra alone. Thus our results can be critical in analyzing 
the implications of present and future y-ray spectral data. 
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2. DETERMINATION OF THE IBL FROM GALAXY 
PHOTON EMISSIVITIES 

Galaxy LFs, <J> V (L), can be determined by properly sampling 
galaxies in a survey for a given wavelength and accounting 
for biases. As in SMS, in this work, we chose to use only 
those references that give values for the integrated LF, i.e., 
the luminosity density (LD), Pl v (z). This is because there is 
generally a lack of knowledge about the covariance of the errors 
in the Schechter function parameters used to determine the 
dominant statistical errors in their analyses. Sufficient numbers 
of galaxies need to be observed up and down the LF for a reliable 
constraint on their total LD, i.e., the integral of the LF. It is not 
even guaranteed that a Schechter function will always provide 
the best analytic fit to the LF data. The co-moving radiation 
energy density, u v (z), is derived from the equivalent co-moving 
specific emissivity £ v (z) = p/.„(z). 

2.1. Luminosity Densities 

In this work, we use the same set of galaxy survey data and 
resulting luminosity densities that covered rest-frame wave- 
lengths from the FUV to the I band as in SMS; however, we 
now extend our calculations to the J and K bands. We have ex- 
cluded the H band due to a paucity of observational data but 
we have checked that our approach is consistent with what data 
exists in this band. We further use the colors derived in the next 
subsection to extend the calculation into the L and M bands (out 
to 4.8 pm). 

The co-moving radiation energy density u v (z) is the time 
integral of the co-moving specific emissivity £„(z), 

/ Zmax 

dz' £Az')—(z'), (1) 

where v' = p(l + z')/(l + z) and z max is the redshift correspond- 
ing to initial galaxy formation (Salamon & Stecker 1998), and 

^(z) = [H 0 (l + zWQ. a + O m (l + z) 3 ] -1 , (2) 

dz 

with £2 a = 0.72 and Ll m = 0.28. We note that we have scaled 
all of the observational data for £ v {z) to a value of h = 0.7 for 
consistency. 


2.2. Average Colors 

The continuum emission from galaxies between 0.8 pm and 
5 pm arises predominantly from stellar photospheres. At these 
near-infrared (NIR) wavelengths, the light is mostly emitted by 
red giant stars. Within a few million years of the first generation 
of star formation, massive stars will have left the main sequence 
to begin populating the red giant branch. At later times, red giant 
branch stars do not vary greatly from one stellar population to 
another. The resulting spectral energy distributions (SEDs) of 
galaxies are therefore much more similar to each other in the 
NIR range than they are at shorter wavelengths. The scatter 
in red and NIR SEDs of galaxies is especially small, since 
all of them contain old, red giant branch stellar populations, 
which dominate the stellar mass and the continuum emission 
around 1 pm. 

Dai et al. (2009) determine galaxy LFs at 3.6, 4.5, 5.8, 
and 8.0 pm for a sample of 3800-5800 galaxies utilizing 
combined photometry from the Spitzer/IRAC Shallow Survey 
with redshifts from the active galactic nucleus and Galaxy 
Evolution Survey of the NOAO Deep Wide-Field Survey Bootes 


Table 1 

Average Colors 


Waveband 

Relative vL v 

Std. Deviation 

Min. Value 

Max. Value 

J (1.2 fim) 

i a 

- 

- 

- 

H (1.6 /xm) 

1.19 

0.08 

1.09 

1.37 

K (2.2 /xm) 

0.97 

0.12 

0.83 

1.23 

L (3.5 /xm) 

0.70 

0.15 

0.43 

0.99 

M (4.8 fim) 

0.50 

0.16 

0.3 

0.84 


Notes. The bias in the colors toward bluer averages occurs because there are not 
many galaxies in the very red categories, so the mean color stays on the blue 
side of (Max + Min)/2. 
a All differences are relative to the J band. 


field. They obtained well defined complete LFs in the local 
redshift bin of z ^ 0.2. They then derived galaxy luminosity 
densities and generated a best-fit SED (their Figure 13). We take 
these data to be representative of low-redshift galaxies (z ^ 0.5). 

To estimate the red and NIR SEDs of galaxies at higher 
redshifts, we have used the template SEDs derived by Kriek 
et al. (2010, 201 1). 4 Those authors utilized extensive multi- 
band photometry obtained for a sample of 3500 K-band selected 
galaxies, at redshifts between z = 0.5 and z = 2.0, close to a 
mass-limited sample of galaxies. For redshifts between z = 0.5 
and z = 1.0, the galaxy stellar masses are uniformly distributed 
between the range of 10 97 and 10 113 M e . At redshifts from 
1 to 2, the mass range covered by these SEDs shrinks to 
10 10 - 3 — 10 1 1 3 . Kriek et al. (2010) grouped their galaxy SEDs into 
30 average templates spanning the full range from the bluest to 
the reddest galaxies. These typical L* galaxies are the ones that 
produce most of the cosmic red and NIR emission. Therefore 
we take the unweighted average of the Kriek et al. (2010, 2011) 
SEDs to represent the full sample of galaxies. This should be 
reasonably accurate as there is only a 0.5 mag Y — L color 
variation (from 1.05 to 3.5 /xm) between the bluest and reddest 
galaxies. Therefore, over this full wavelength range, the average 
colors we have adopted must be quite close to the colors of any 
galaxy, at least at these redshifts. We express the average colors 
we have adopted as flux ratios normalized to J band in Table 1 . 

The peak in the galaxy SEDs that we used occurs around a rest 
wavelength of 1.6 pm. We note that there is a subtle trend for 
the red and NIR SEDs to become bluer at higher redshift. This 
is because the younger stellar populations were more important 
in the early universe. This effect produces a small shift of the 
SED peak to slightly shorter rest wavelengths at higher redshifts. 
This trend may continue to redshifts higher than two. However, 
because this effect is small, we assume that the SEDs of Kriek 
et al. (2010) apply to all higher redshifts. We assume that these 
average colors also apply to the less massive galaxies that Kriek 
et al. could not study. Fortunately, the total NIR emission of 
galaxies tends to be dominated by the more massive galaxies 
(e.g., Ly et al. 2011). Furthermore, galaxies of masses smaller 
by an order of magnitude are only slightly bluer. 


2.3. Photon Density Calculations 

As in SMS, we derive a luminosity confidence band in each 
of our additional wavebands by using a robust rational fitting 
function characterized by 


Pl„ = £ v (z) = 


ax + b 
cx 2 + dx + e 


(3) 


4 See also http://astro.berkeley.edu/~mariska/comp/. 
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Figure 1. Specific emissivities for J, K, L, andM wavebands. In the L and 
text in order to fully determine the specific emissivities in these wavebands. 


M panels, the ./- and A'-band data have been shifted using the color relations given in the 
The gray shading represents the 68% confidence bands (see text). 


where x = log(l + z) and a, b, c, d, and e are free parameters. 
We compute the 68% confidence band from Monte Carlo 
simulations by finding 100,000 realizations of the data and then 
fitting the rational function given by Equation (3). As in SMS, 
we model symmetric error bars with a Gaussian distribution, 
while choosing a skew normal distribution to model asymmetric 
errors. 

The A'-band data do not extend beyond a redshift of z ~ 2, 
while the ./-band data extend to z ~ 3. Thus, the fits determined 
using the above method can not be trusted beyond these limits. 
To handle this issue, we truncate our fits at a redshift of z — 2, 
and then color transform our previous I band from SMS to fit at 
the higher redshifts. This is a reasonable assumption owing 
to the expected similarity in the SEDs of galaxies at these 
wavelengths (see previous sections). Indeed, the overlapping 
7-band data beyond a redshift of two are in good agreement 
with this color-shifted band. 

Figure 1 shows the redshift evolution of the LD for the various 
wavebands based on the survey data published in the literature 
for the wavelength bands that we use to extend the previous work 
of SMS. 5 The upper and lower limits of the bands correspond 
to the highest and lowest IBL consistent with 68% confidence 
errors. The error bars on the 7- and A'-band data taken from 
Pozzetti et al. (2003) are small in comparison to other data 
used. Their uncertainty was estimated strictly by considering the 
range of acceptable Schechter parameters when fitting their LFs. 
As other authors considered systematic errors such as cosmic 
variance and extrapolation to the faint end, we have chosen to 
assign a linear error of 20% to these data to bring their errors 

5 References for the values of £ v (z) used to construct Figure 1 are listed in 
the legend. 


more in line with what is typically quoted so that these data 
do not inordinately influence the fitting routine. The upper and 
lower limits of the bands correspond to the highest and lowest 
IBL consistent with 68% confidence errors. 

Our method of dealing with the confidence band for redshifts 
beyond where there are any data for 7 and K is not critical for the 
opacity calculation. There is very little contribution to the y-ray 
opacities from photons in these wavebands beyond a redshift of 
2 because of the short time interval of the emission from galaxies 
at higher redshifts (see Equation (2)). To verify this, we set an 
upper limit by assuming that the value of the photon density at 
z = 2 remains constant out to higher redshifts. Since we expect 
that the LD will drop off at the higher redshifts, this assumption 
gives an upper limit. We found no appreciable differences in the 
calculated opacities under this assumption (see next section) as 
expected. 

Having determined the colors between the 7, K , L, and 
M bands, we use them to transform the 7 and K LDs into the L 
and M bands and again apply our fitting technique to produce 
confidence limits for those wavelengths. The color- transformed 
data along with the confidence bands for L and M can also be 
seen in Figure 1 . We note that the wider resulting confidence 
bands for L and M occur because we have propagated the error 
determined from the colors to the 7 and K data (reflected in 
the size of the error bars) before running the Monte Carlo 
simulations and fits. This is as expected, since the band of 
acceptable photon densities derived from the data should reflect 
both errors. As in SMS, we then interpolate between the upper 
limit and lower limit of the confidence bands between the various 
wavebands separately to find the upper and lower limit rest- 
frame luminosity densities. 
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Figure 2. Empirically determined opacities for redshifts of 0.1, 0.5, 1, 3, and 5 
extended from SMS. The dashed lines indicate the opacities r = 1 and r = 3. 


3. DETERMINATION OF THE OPTICAL DEPTH TO 
y-RAYS AND RESULTING y-RAY HORIZON 

Using the results of the previous section, the specific emis- 
sivity can then be derived for the highest and lowest IBL LDs. 
The co-moving radiation energy density is then determined from 
Equation (1). The photon densities derived thereby are given by 


n(e, z) = u(e, z)/e 


(4) 


with e = hv, and with u(e, z ) given by Equation (1). 

With the co-moving photon density n(e, z) evaluated, the opti- 
cal depth for y- rays owing to electron-positron pair-production 
interactions with photons of the stellar radiation background can 
be determined from the expression (Stecker et al. 1992) 

r z ' dt r 2 x 

x(Eq, z e )=c dz — dx - 
Jo dz Jo 2 

fOO 

x / de(\+zfn{e,z)oyy[s{z)]. (5) 
Jtth 

In Equation (5), dt/dz is given by Equation (2), Eq is the 
observed y- ray energy at redshift zero, e is evaluated at redshift 
Z, z e is the redshift of the y-ray source, x = (1 — cos (9), 0 
being the angle between the y-ray and the soft background 
photon, and the pair-production cross section a YY is zero for the 
rest system center-of-mass energy ^fs < 2m e c 2 , m e being the 
electron mass. Above this threshold, the pair-production cross 
section is given by 


o- rr (s)= T^cr-rC 1 - ft~) 


2 ft (ft 2 — 2) + (3 — ft 4 ) In 


1 + i 


( 6 ) 


where oj is the Thompson scattering cross section and ft = 
(1 — 4/71 2 c 4 /s) 1 ^ 2 (Breit & Wheeler 1934; Jauch & Rohrlich 
1955). 



Figure 3. A r = 1 energy-redshift plot (Fazio & Stecker 1970) showing our 
uncertainty band results compared with the Fermi plot of their highest energy 
photons from FSRQs (red), BL Lacs (black), and y-ray bursts (blue) vs. redshift 
(from Abdo et al. 2010). 

(A color version of this figure is available in the online journal.) 


As derived in SMS, the pair-production cross section energy 
has a threshold at k = 4.8 //m • £ y (TeV), determined from the 
energy required to produce twice the electron rest mass in the 
center of mass frame. Since the maximum X is in the rest-frame 
Mband at 4.8 // m at redshift z, the resulting energy at threshold 
is ~ 1 .6 Te V at’ = 0. The energy at interaction in the rest frame is 
given by (1 +z)E y meaning the maximum y-ray energy affected 
by the photon range that we consider is ~1 .6 (1 + z) _1 TeV. 

The 68% opacity ranges for z = 0.1, 0.5, 1,3 and 5, 
calculated using the SMS methods as described above, are 
plotted in Figure 2. The increasing uncertainties in the y-ray 
opacity toward higher redshifts are a reflection of the increasing 
widths of the uncertainty bands in the luminosity densities 
shown in Figure 1, as follows from Equations (4) and (5). 

Figure 3 shows an energy-redshift plot of the highest energy 
photons from extragalactic sources at various redshifts from 
Fermi as given by Abdo et al. (2010) along with our 68% 
confidence band for r = 1, extending our result from SMS 
down to a redshift of z = 0.2. We stress that our 68% confidence 
band as shown in Figure 3 is the range of optical depths allowed 
as derived from our observationally determined IBL results. 

4. APPLICATION OF OUR RESULTS TO PKS 1424+240 

To show the utility of our results, we consider redshift 
limits obtained from y-ray observations of the distant blazar 
PKS 1424+240. BL Lac objects such as this one typically dis- 
play bright and featureless continuous spectra making definite 
spectroscopic determinations of their redshifts challenging. A 
recent determination of a lower limit of the redshift of PKS 
1424+240 has been obtained by Furniss et al. (2013). They find 
z Js 0.6035, inferred from Hubble observations of the measured 
positions of Lyman lines from intervening Lya absorbers at 
lower redshifts. PKS 1424+240 has been observed at y-ray en- 
ergies lower than 100 GeV by Fermi - LAT and at energies above 
500 GeV by VERITAS (Acciari et al. 2010). 

Using our opacity results, we can place an upper limit on the 
redshift of this source. A power-law spectrum with index ~1.8 
has been determined from Fermi for energies below 100 GeV 
(Furniss et al. 2013). In this energy range, there should be 
little y-ray attenuation. Therefore, extending this power-law 
spectrum out to energies measured by VERITAS should represent 
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Figure 4. Determination of the upper limit redshift for PKS 1424+240. The 
points represent the upper limits to the y - ray opacity derived from the VERITAS 
data. The dashed curve corresponds to our lower limit opacity for z = .6035 
while the solid curve is our best-fit lower limit opacity to the opacity limits 
derived from the results shown in Figure 3, corresponding to an upper limit 
redshift of z = 1.0. 

an upper limit on its intrinsic spectrum. An upper limit for 
the opacity for each of the measured /-ray energies can then 
be determined from the flux difference between the extended 
power-law spectrum and the measured fluxes. This upper limit 
can then be compared with our redshift-dependent opacity 
curves generated from our lower limit IBL, thereby determining 
an upper limit on the redshift of PKS 1424+240. Our upper limit 
on z is found to be 1 .0; our lower limit opacities for redshifts 
greater than z = 1.0 exceed the upper limits on the opacity 
determined from the flux differences. Figure 4 shows the upper 
limits on the opacities determined from the VERITAS data along 
with our lower limit z = 1 0 opacity curve. Aleksic et al. (2014) 
have used MAGIC data to place an upper limit on the redshift 
of PKS 1424+240 of 0.81. However, the limit obtained by the 
MAGIC collaboration depends on assuming absolute accuracy 
of the IBL model of Franceschini et al. (2008). The loosening of 
this limit to a slightly larger upper limit of z — 1 .0 is a result of 
using our observationally derived uncertainty band rather than 
a theoretical model. 

5. COMPARISON WITH OTHER RESULTS 

Helgason & Kashlinsky (2012, HK) have used a method sim- 
ilar to ours, using galaxy LFs to determine the /-ray opacity 
of the universe with uncertainties over the same energy range 
as considered here. Their results are also model independent. 
However, their treatment differs from our ours in two respects: 
( 1 ) the HK uncertainty band is determined by fitting to Schechter 
function parameterizations with upper and lower limits deter- 
mined by fitting to faint end slopes of the Schechter functions; 
we use observationally determined LDs with observational er- 
rors combined with our Monte Carlo treatment as described in 
Section 2.3. Our use of observationally determined errors takes 
account of systematic effects such as cosmic variance. (2) We 
use color data as described in Section 2.3 to extrapolate our LDs 
to higher redshifts, allowing us to make use of more observa- 
tional data, particularly in the L and M bands; HK postulate an 
exponential cutoff in their LDs at higher redshifts. The result of 
these differences is that our opacities are somewhat larger than 
those of HK. 

Our uncertainty band for the IBL at z — 0 (the EBL) is in 
full agreement with that obtained by Abramowski et al. (2013) 
using /-ray spectra from H.E.S.S. (see their Figure 5). Our 


results are consistent with the work of the Fermi collaboration 
(2010) as shown in our Figure 3. Our results for z = 1 are in fair 
agreement with those of the analyses of the Fermi collaboration 
(Ackermann et al. 2012). Our r = 1 (horizon) band is consistent 
with the result of Dominguez et al. (2013). 

6. CONCLUSIONS 

We have presented an extension of our previous determination 
of 68% confidence bands giving upper and lower limits on the 
IBL out to 4.8 fim. This model-independent determination is 
entirely based on observationally derived LFs from local and 
deep galaxy survey data and color data. This has enabled us 
to directly derive both the /-ray opacity and its observational 
uncertainties as a function of both energy and redshift out to an 
energy ofl.6/(l+z) Te V for z ^ 5 . We have applied our lower 
limit opacities to the results of both ground-based air Cerenkov 
telescope and Fermi-LKV observations of PKS 1424+240, 
allowing us to place a new upper limit on the redshift of this 
source, viz., z ^ 1.0, independent of IBL modeling. 

We find no direct evidence in our spectral analysis for a re- 
quired modification of our predicted / -ray opacities as shown in 
Figure 4, and for the opacity in the direction of PKS 1424+240, 
either owing to axion-photon mixing as suggested Meyer & 
Horns (2013) or by the existence of secondary production 
effects as suggested by Essey & Kusenko (2013). However, 
such effects may need to be invoked should a subsequently 
determined redshift of PKS 1424+240 turn out to be greater 
than 1. 

We thank Amy Furniss for supplying us with the VERITAS 
/-ray data on PKS 1424+240 that we used in our analysis. We 
also thank Mariska Kriek for supplying us with galaxy template 
SEDs. 
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